"""
    RadianceBUFRReader

A Julia module for reading satellite radiance BUFR files used in GSI.

This module provides functionality for reading satellite brightness temperature
observations from BUFR format files (e.g., gdas1.t12z.1bamua.tm00.bufr_d).

Unlike PrepBUFR (conventional observations), radiance BUFR files contain:
- Satellite brightness temperatures (observed radiances)
- Multi-channel observations (e.g., 15 channels for AMSU-A)
- Satellite geometry (zenith/azimuth angles)
- Sensor metadata (satellite ID, instrument type)

# Supported Instruments

## Microwave Sounders
- **AMSU-A** (Advanced Microwave Sounding Unit-A): 15 channels, NOAA-15-19, MetOp
- **MHS** (Microwave Humidity Sounder): 5 channels, NOAA-18-19, MetOp
- **ATMS** (Advanced Technology Microwave Sounder): 22 channels, SNPP, NOAA-20

## Infrared Sounders
- **IASI** (Infrared Atmospheric Sounding Interferometer): 616 channels, MetOp
- **AIRS** (Atmospheric Infrared Sounder): 2378 channels, Aqua
- **CrIS** (Cross-track Infrared Sounder): 1305 channels, SNPP, NOAA-20

## Imagers
- **HIRS** (High-resolution Infrared Radiation Sounder): 20 channels
- **AVHRR** (Advanced Very High Resolution Radiometer): 6 channels

# Usage

```julia
using GSICoreAnalysis.DataIO.RadianceBUFRReader

# Read AMSU-A radiance observations
radiance_obs = read_radiance_bufr("/path/to/gdas1.t12z.1bamua.tm00.bufr_d")

# Access data
println("Satellite: ", radiance_obs.satellite_id)
println("Instrument: ", radiance_obs.instrument)
println("Channels: ", radiance_obs.channels)
println("Observations: ", size(radiance_obs.brightness_temp))
println("Locations: ", length(radiance_obs.latitude))
```

# Data Structure

The `RadianceObservation` structure contains:
- `satellite_id`: Satellite identifier (e.g., "noaa-15")
- `instrument`: Instrument name (e.g., "amsua")
- `channels`: Vector of channel numbers
- `brightness_temp`: Matrix of brightness temperatures [K] (n_channels × n_obs)
- `latitude`, `longitude`: Observation geolocation
- `zenith_angle`, `azimuth_angle`: Satellite viewing geometry
- `solar_zenith`, `solar_azimuth`: Solar geometry
- `qc_flags`: Quality control indicators
- `obs_time`: Observation times

# Integration with CRTM

The output format is designed to integrate seamlessly with CRTM forward/adjoint operators:

```julia
# Read observations
obs = read_radiance_bufr("gdas1.t12z.1bamua.tm00.bufr_d")

# Setup CRTM
crtm_init_julia("amsua_n15", coef_path, obs.channels)

# Forward operator
tb_simulated = crtm_forward_julia(temp, pres, humid, skin_temp, wind_speed,
                                   obs.zenith_angle, obs.azimuth_angle, obs.channels)

# Innovation
innovation = obs.brightness_temp .- tb_simulated
```

# References
- NCEP BUFR Tables: https://www.emc.ncep.noaa.gov/mmb/data_processing/NCEP_BUFR_File_Structure.htm
- GSI radiance processing: comGSIv3.7/src/read_radiag.f90
"""
module RadianceBUFRReader

using Dates
using Printf
using Statistics
using eccodes_jll

using ..RadianceInstrumentConfig: instrument_info

# Import parent types
using ..DataIO: ObservationData

export RadianceObservation, read_radiance_bufr, RadianceBUFRConfig
export get_satellite_name, get_instrument_name

"""
    RadianceObservation

Structure containing satellite radiance observation data.

# Fields
- `satellite_id::String`: Satellite identifier (e.g., "noaa-15", "metop-a")
- `instrument::String`: Instrument name (e.g., "amsua", "mhs", "iasi")
- `channels::Vector{Int32}`: Channel numbers
- `brightness_temp::Matrix{Float32}`: Brightness temperatures [K] (n_channels × n_obs)
- `latitude::Vector{Float32}`: Observation latitudes [degrees]
- `longitude::Vector{Float32}`: Observation longitudes [degrees]
- `zenith_angle::Vector{Float32}`: Satellite zenith angle [degrees]
- `azimuth_angle::Vector{Float32}`: Satellite azimuth angle [degrees]
- `solar_zenith::Vector{Float32}`: Solar zenith angle [degrees]
- `solar_azimuth::Vector{Float32}`: Solar azimuth angle [degrees]
- `qc_flags::Matrix{Int32}`: Quality control flags (n_channels × n_obs)
- `obs_error::Matrix{Float32}`: Observation errors [K] (n_channels × n_obs)
- `obs_time::Vector{DateTime}`: Observation times
- `metadata::Dict{String,Any}`: Additional metadata
"""
struct RadianceObservation
    satellite_id::String
    instrument::String
    channels::Vector{Int32}
    brightness_temp::Matrix{Float32}
    latitude::Vector{Float32}
    longitude::Vector{Float32}
    zenith_angle::Vector{Float32}
    azimuth_angle::Vector{Float32}
    solar_zenith::Vector{Float32}
    solar_azimuth::Vector{Float32}
    qc_flags::Matrix{Int32}
    obs_error::Matrix{Float32}
    obs_time::Vector{DateTime}
    metadata::Dict{String,Any}
end

"""
    RadianceBUFRConfig

Configuration for radiance BUFR reading.

# Fields
- `max_obs::Union{Int,Nothing}`: Maximum observations to read (nothing = all)
- `channel_subset::Union{Vector{Int},Nothing}`: Specific channels to extract (nothing = all)
- `qc_threshold::Int`: Maximum QC flag value to accept
- `angle_limit::Float64`: Maximum satellite zenith angle [degrees]
- `verbose::Bool`: Print diagnostic information
"""
struct RadianceBUFRConfig
    max_obs::Union{Int,Nothing}
    channel_subset::Union{Vector{Int},Nothing}
    qc_threshold::Int
    angle_limit::Float64
    verbose::Bool

    function RadianceBUFRConfig(;
        max_obs::Union{Int,Nothing}=nothing,
        channel_subset::Union{Vector{Int},Nothing}=nothing,
        qc_threshold::Int=3,
        angle_limit::Float64=65.0,
        verbose::Bool=false
    )
        new(max_obs, channel_subset, qc_threshold, angle_limit, verbose)
    end
end

# Satellite ID lookup table (WMO code → name)
const SATELLITE_IDS = Dict(
    1 => "noaa-15",
    2 => "noaa-16",
    3 => "noaa-17",
    4 => "noaa-18",
    5 => "noaa-19",
    206 => "metop-a",
    207 => "metop-b",
    208 => "metop-c",
    224 => "snpp",
    225 => "noaa-20",
    226 => "noaa-21",
    783 => "aqua"
)

# Instrument ID lookup table (WMO code → name)
const INSTRUMENT_IDS = Dict(
    570 => "amsua",      # AMSU-A
    574 => "mhs",        # MHS
    203 => "mhs",        # MHS (NOAA local code)
    420 => "atms",       # ATMS
    221 => "iasi",       # IASI
    591 => "airs",       # AIRS
    620 => "cris",       # CrIS
    606 => "hirs",       # HIRS
    207 => "avhrr"       # AVHRR
)

# Default observation errors by instrument [K]
const DEFAULT_OBS_ERRORS = Dict(
    "amsua" => collect(instrument_info("amsua").error_std),
    "mhs" => collect(instrument_info("mhs").error_std),
    "atms" => collect(instrument_info("atms").error_std),
    "iasi" => fill(1.5, 616),
    "airs" => fill(2.0, 2378),
    "cris" => fill(1.5, 1305),
    "hirs" => fill(2.0, 20),
    "avhrr" => fill(2.0, 6)
)

"""
    get_satellite_name(said::Int) -> String

Get satellite name from WMO satellite identifier.
"""
function get_satellite_name(said::Int)
    return get(SATELLITE_IDS, said, "unknown-$said")
end

"""
    get_instrument_name(siid::Int) -> String

Get instrument name from WMO instrument identifier.
"""
function get_instrument_name(siid::Int)
    return get(INSTRUMENT_IDS, siid, "unknown-$siid")
end

"""
    read_radiance_bufr(filename::String; config::RadianceBUFRConfig=RadianceBUFRConfig()) -> RadianceObservation

Read satellite radiance observations from BUFR file.

# Arguments
- `filename`: Path to radiance BUFR file (e.g., gdas1.t12z.1bamua.tm00.bufr_d)
- `config`: Configuration options

# Returns
- `RadianceObservation`: Structure containing brightness temperatures and metadata

# Example
```julia
obs = read_radiance_bufr("/path/to/gdas1.t12z.1bamua.tm00.bufr_d")
println("Read \$(size(obs.brightness_temp, 2)) observations with \$(length(obs.channels)) channels")
```
"""
function read_radiance_bufr(filename::String; config::RadianceBUFRConfig=RadianceBUFRConfig())
    if !isfile(filename)
        error("Radiance BUFR file not found: $filename")
    end

    config.verbose && @info "Reading radiance BUFR file: $filename"

    return read_radiance_bufr_eccodes(filename, config)
end

"""
    read_radiance_bufr_eccodes(filename::String, config::RadianceBUFRConfig) -> RadianceObservation

Read radiance BUFR files using the native ecCodes bindings.
"""
const ECCODES_LIB = eccodes_jll.eccodes
const ECCODES_MISSING_THRESHOLD = -1.0e99

const PLANCK_C1 = 1.1910427e-8  # W·cm²·sr⁻¹
const PLANCK_C2 = 1.4387752     # cm·K

function read_radiance_bufr_eccodes(filename::String, config::RadianceBUFRConfig)
    _configure_eccodes_environment!()

    observations = Vector{Dict{String,Any}}()
    file = ccall(:fopen, Ptr{Cvoid}, (Cstring, Cstring), filename, "rb")
    file == C_NULL && error("Failed to open radiance BUFR file: $filename")
    max_obs = config.max_obs

    try
        while true
            err = Ref{Cint}(0)
            handle = ccall((:codes_bufr_handle_new_from_file, ECCODES_LIB), Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{Cvoid}, Ref{Cint}), C_NULL, file, err)
            handle == C_NULL && break

            try
                category = _get_scalar_long(handle, "dataCategory"; default=-1)
                ccall((:codes_set_long, ECCODES_LIB), Cint, (Ptr{Cvoid}, Cstring, Clong), handle, "unpack", 1)
                subsets = _get_scalar_long(handle, "numberOfSubsets"; default=0)

                if category == 21 && subsets > 0
                    _decode_radiance_message!(handle, config, observations, max_obs)
                    if max_obs !== nothing && length(observations) >= max_obs
                        break
                    end
                end
            finally
                ccall((:codes_handle_delete, ECCODES_LIB), Cint, (Ptr{Cvoid},), handle)
            end
        end
    finally
        ccall(:fclose, Cint, (Ptr{Cvoid},), file)
    end

    if isempty(observations)
        error("No radiance observations found in file: $filename")
    end

    if max_obs !== nothing && length(observations) > max_obs
        resize!(observations, max_obs)
    end

    config.verbose && @info "Decoded $(length(observations)) observation locations (ecCodes)"

    return convert_radiance_observations(observations, filename, config)
end

struct ChannelData
    channels::Vector{Int}
    brightness::Vector{Vector{Float64}}
    corrections::Vector{Float64}
    source_kind::String
end

function _decode_radiance_message!(handle::Ptr{Cvoid}, config::RadianceBUFRConfig, observations::Vector{Dict{String,Any}}, max_obs::Union{Nothing,Int})
    lat_candidates = (
        "#1#latitude",
        "#1#latitudeHighAccuracy",
        "#1#latitudeCoarseAccuracy",
        "#1#geographicalLatitude",
        "#1#latitudeOfObservation",
    )
    lon_candidates = (
        "#1#longitude",
        "#1#longitudeHighAccuracy",
        "#1#longitudeCoarseAccuracy",
        "#1#geographicalLongitude",
        "#1#longitudeOfObservation",
    )

    lat, lat_key = _select_vector(handle, lat_candidates; default=NaN)
    n_points = length(lat)
    if n_points == 0 || all(isnan, lat)
        config.verbose && @warn "Radiance message missing latitude; skipping BUFR subset" lat_key=lat_key
        return
    end
    lon, lon_key = _select_vector(handle, lon_candidates; expected=n_points, default=NaN)
    if isempty(lon) || all(isnan, lon)
        config.verbose && @warn "Radiance message missing longitude; skipping BUFR subset" lat_key=lat_key lon_key=lon_key
        return
    end
    # instrument-specific diagnostics can be added here if needed

    zenith = _get_vector(handle, "#1#satelliteZenithAngle"; expected=n_points, default=0.0)
    azimuth = _get_vector(handle, "#1#bearingOrAzimuth"; expected=n_points, default=0.0)
    solar_zen = _get_vector(handle, "#1#solarZenithAngle"; expected=n_points, default=0.0)
    solar_az = _get_vector(handle, "#1#solarAzimuth"; expected=n_points, default=0.0)
    minute = _get_vector(handle, "#1#minute"; expected=n_points, default=0.0)
    second = _get_vector(handle, "#1#second"; expected=n_points, default=0.0)

    year = _get_scalar_long(handle, "#1#year"; default=2000)
    month = _get_scalar_long(handle, "#1#month"; default=1)
    day = _get_scalar_long(handle, "#1#day"; default=1)
    hour = _get_scalar_long(handle, "#1#hour"; default=0)

    sat_id = _get_scalar_long(handle, "#1#satelliteIdentifier"; default=0)
    inst_id = _get_scalar_long(handle, "#1#satelliteInstruments"; default=0)
    instrument_name = get_instrument_name(inst_id)

    channel_data = _collect_channel_data(handle, n_points, instrument_name)
    channel_ids = channel_data.channels
    if isempty(channel_ids)
        return
    end
    bt_arrays = channel_data.brightness

    limit = max_obs === nothing ? n_points : min(n_points, max_obs - length(observations))
    limit <= 0 && return

    for idx in 1:limit
        brightness = Vector{Float64}(undef, length(channel_ids))
        for (ch_idx, arr) in enumerate(bt_arrays)
            brightness[ch_idx] = _sanitize_bt(arr[idx])
        end

        obs = Dict{String,Any}(
            "latitude" => lat[idx],
            "longitude" => lon[idx],
            "zenith_angle" => zenith[idx],
            "azimuth_angle" => azimuth[idx],
            "solar_zenith" => solar_zen[idx],
            "solar_azimuth" => solar_az[idx],
            "channels" => channel_ids,
            "brightness_temps" => brightness,
            "satellite_id" => sat_id,
            "instrument_id" => inst_id,
            "source_value_kind" => channel_data.source_kind,
            "base_time" => (
                year,
                month,
                day,
                hour,
                Int(clamp(round(minute[idx]), 0, 59)),
                Int(clamp(round(second[idx]), 0, 59)),
            ),
        )

        push!(observations, obs)
        if max_obs !== nothing && length(observations) >= max_obs
            break
        end
    end
end

function _select_vector(handle::Ptr{Cvoid}, candidates::Tuple{Vararg{String}}; expected::Union{Nothing,Int}=nothing, default::Float64=NaN)
    for key in candidates
        vec = _get_vector(handle, key; expected=expected, default=default, required=false)
        isempty(vec) && continue
        if any(!isnan, vec)
            return vec, key
        end
    end
    return Float64[], nothing
end

function _configure_eccodes_environment!()
    overlay = normpath(joinpath(@__DIR__, "..", "..", "deps", "eccodes_overrides"))
    paths = String[]
    if isdir(overlay)
        push!(paths, overlay)
    end

    eccodes_home = get(ENV, "ECCODES_HOME", nothing)
    if eccodes_home !== nothing
        defs_path = joinpath(eccodes_home, "share", "eccodes", "definitions")
        if isdir(defs_path)
            push!(paths, defs_path)
        end
    end

    libdir = dirname(ECCODES_LIB)
    artifact_defs = joinpath(dirname(libdir), "share", "eccodes", "definitions")
    if isdir(artifact_defs)
        push!(paths, artifact_defs)
    end

    existing = get(ENV, "ECCODES_DEFINITION_PATH", "")
    for part in split(existing, ':')
        if !isempty(part)
            push!(paths, part)
        end
    end

    unique_defs = unique(filter(isdir, paths))
    if !isempty(unique_defs)
        ENV["ECCODES_DEFINITION_PATH"] = join(unique_defs, ":")
    end

    sample_paths = String[]
    if eccodes_home !== nothing
        samples = joinpath(eccodes_home, "share", "eccodes", "samples")
        if isdir(samples)
            push!(sample_paths, samples)
        end
    end

    artifact_samples = joinpath(dirname(libdir), "share", "eccodes", "samples")
    if isdir(artifact_samples)
        push!(sample_paths, artifact_samples)
    end

    existing_samples = get(ENV, "ECCODES_SAMPLES_PATH", "")
    for part in split(existing_samples, ':')
        if !isempty(part)
            push!(sample_paths, part)
        end
    end

    unique_samples = unique(filter(isdir, sample_paths))
    if !isempty(unique_samples)
        ENV["ECCODES_SAMPLES_PATH"] = join(unique_samples, ":")
    end
end

function _collect_channel_data(handle::Ptr{Cvoid}, n_points::Int, instrument::String)
    channels = Int[]
    brightness = Vector{Vector{Float64}}()
    corrections = Float64[]
    subset = 1
    used_radiance = false


    while true
        ch_key = "#" * string(subset) * "#channelNumber"
        size_ref = Ref{Clong}(0)
        rc = ccall((:codes_get_size, ECCODES_LIB), Cint, (Ptr{Cvoid}, Cstring, Ref{Clong}), handle, ch_key, size_ref)
        rc != 0 && break

        vals = Vector{Clong}(undef, Int(size_ref[]))
        ccall((:codes_get_long_array, ECCODES_LIB), Cint, (Ptr{Cvoid}, Cstring, Ptr{Clong}, Ref{Clong}), handle, ch_key, pointer(vals), size_ref)
        push!(channels, Int(vals[1]))

        bt_key = "#" * string(subset) * "#brightnessTemperature"
        bt = _get_vector(handle, bt_key; expected=n_points, default=NaN, fill_nan=true)

        if all(x -> isnan(x), bt)
            rad_key = "#" * string(subset) * "#channelRadiance"
            radiance = _get_vector(handle, rad_key; expected=n_points, default=NaN, fill_nan=true)
            if any(isfinite, radiance)
                used_radiance = true
                wn = _instrument_wavenumber(instrument, Int(vals[1]))
                bt = [_radiance_to_brightness(instrument, wn, r) for r in radiance]
            end
        end

        push!(brightness, bt)

        corr_key = "#" * string(subset) * "#coldSpaceTempCorrection"
        corr_vec = _get_vector(handle, corr_key; expected=1, default=NaN)
        push!(corrections, corr_vec[1])

        subset += 1
    end

    source_kind = used_radiance ? "channelRadiance" : "brightnessTemperature"
    return ChannelData(channels, brightness, corrections, source_kind)
end

function _instrument_wavenumber(instrument::AbstractString, channel::Int)
    if instrument == "cris"
        return _cris_wavenumber(channel)
    else
        return NaN
    end
end

function _cris_wavenumber(channel::Int)
    if channel <= 0
        return NaN
    elseif channel <= 713
        return 650.0 + 0.625 * (channel - 1)
    elseif channel >= 714 && channel <= 1570
        return 1210.0 + 0.625 * (channel - 714)
    elseif channel >= 1571
        return 2165.625 + 0.625 * (channel - 1596)
    else
        return NaN
    end
end

function _radiance_to_brightness(instrument::AbstractString, wavenumber_cm::Float64, radiance::Float64)
    if !isfinite(radiance) || radiance <= 0.0
        return NaN
    end
    if instrument == "cris"
        if !isfinite(wavenumber_cm) || wavenumber_cm <= 0.0
            return NaN
        end
        ratio = PLANCK_C1 * wavenumber_cm^3 / radiance
        if !isfinite(ratio) || ratio <= 0.0
            return NaN
        end
        denom = log1p(ratio)
        if !isfinite(denom) || denom <= 0.0
            return NaN
        end
        return PLANCK_C2 * wavenumber_cm / denom
    end
    return NaN
end

function _get_vector(handle::Ptr{Cvoid}, key::String; expected::Union{Int,Nothing}=nothing, default::Float64=NaN, required::Bool=false, fill_nan::Bool=false)
    size_ref = Ref{Clong}(0)
    rc = ccall((:codes_get_size, ECCODES_LIB), Cint, (Ptr{Cvoid}, Cstring, Ref{Clong}), handle, key, size_ref)

    if rc != 0
        if required
            error("Missing required key '$key' in BUFR message")
        end
        len = expected === nothing ? 0 : expected
        return fill(default, len)
    end

    len = Int(size_ref[])
    arr = Vector{Float64}(undef, len)
    ccall((:codes_get_double_array, ECCODES_LIB), Cint, (Ptr{Cvoid}, Cstring, Ptr{Float64}, Ref{Clong}), handle, key, pointer(arr), size_ref)

    target_len = expected === nothing ? len : expected
    if len == target_len
        return fill_nan ? _sanitize_bt_array(arr) : arr
    elseif len == 1 && target_len > 1
        vec = fill(arr[1], target_len)
        return fill_nan ? _sanitize_bt_array(vec) : vec
    elseif len > target_len && target_len > 0
        return fill_nan ? _sanitize_bt_array(arr[1:target_len]) : arr[1:target_len]
    else
        vec = fill(default, target_len)
        copy_len = min(len, target_len)
        if copy_len > 0
            vec[1:copy_len] .= arr[1:copy_len]
        end
        return fill_nan ? _sanitize_bt_array(vec) : vec
    end
end

function _sanitize_bt_array(arr::Vector{Float64})
    for i in eachindex(arr)
        arr[i] = _sanitize_bt(arr[i])
    end
    return arr
end

function _sanitize_bt(val::Float64)
    if !isfinite(val) || val <= ECCODES_MISSING_THRESHOLD
        return NaN
    end
    return val
end

function _get_scalar_long(handle::Ptr{Cvoid}, key::String; default::Int=0)
    value = Ref{Clong}(0)
    rc = ccall((:codes_get_long, ECCODES_LIB), Cint, (Ptr{Cvoid}, Cstring, Ref{Clong}), handle, key, value)
    if rc != 0
        return default
    end
    return Int(value[])
end

"""
    convert_radiance_observations(observations::Vector{Dict}, filename::String, config::RadianceBUFRConfig) -> RadianceObservation

Convert raw observation dictionaries to RadianceObservation structure.
"""
function convert_radiance_observations(observations::Vector{Dict{String,Any}}, filename::String, config::RadianceBUFRConfig)
    nobs = length(observations)

    # Determine instrument and satellite from first observation
    first_obs = observations[1]
    satellite_id_code = get(first_obs, "satellite_id", 0)
    instrument_id_code = get(first_obs, "instrument_id", 0)
    source_value_kind = get(first_obs, "source_value_kind", "brightnessTemperature")

    satellite_name = get_satellite_name(satellite_id_code)
    instrument_name = get_instrument_name(instrument_id_code)

    # Get channel list from first observation
    channels_list = get(first_obs, "channels", Int[])
    nchans = length(channels_list)
    channels = convert(Vector{Int32}, channels_list)

    settings = nothing
    try
        candidate = instrument_info(instrument_name)
        if length(candidate.channels) == nchans
            settings = candidate
        else
            @warn "Channel count mismatch" instrument=instrument_name expected=length(candidate.channels) found=nchans
        end
    catch
        settings = nothing
    end

    config.verbose && @info "Instrument: $instrument_name on $satellite_name, $nchans channels"

    # Apply channel subset filter if specified
    full_channel_count = length(channels_list)

    if config.channel_subset !== nothing
        channel_indices = findall(in(config.channel_subset), channels)
        channels = channels[channel_indices]
        nchans = length(channels)
    else
        channel_indices = 1:nchans
    end

    # Allocate output arrays
    brightness_temp = zeros(Float32, nchans, nobs)
    latitude = zeros(Float32, nobs)
    longitude = zeros(Float32, nobs)
    zenith_angle = zeros(Float32, nobs)
    azimuth_angle = zeros(Float32, nobs)
    solar_zenith = zeros(Float32, nobs)
    solar_azimuth = zeros(Float32, nobs)
    qc_flags = zeros(Int32, nchans, nobs)
    obs_time = Vector{DateTime}(undef, nobs)

    # Get default observation errors and use flags for this instrument
    if settings !== nothing && length(settings.error_std) == full_channel_count
        base_errors = Float32.(settings.error_std)
        use_flags = collect(settings.use_flags)
    else
        fallback = get(DEFAULT_OBS_ERRORS, instrument_name, fill(2.0, full_channel_count))
        base_errors = Float32.(fallback)
        use_flags = fill(Int8(1), full_channel_count)
    end

    if length(base_errors) < full_channel_count
        base_errors = vcat(base_errors, fill(Float32(2.0), full_channel_count - length(base_errors)))
    end
    if length(use_flags) < full_channel_count
        use_flags = vcat(use_flags, fill(Int8(1), full_channel_count - length(use_flags)))
    end

    if settings !== nothing && length(settings.tb_min) == full_channel_count
        base_tb_min = Float32.(settings.tb_min)
        base_tb_max = Float32.(settings.tb_max)
        required_valid_fraction = settings.min_valid_fraction
    else
        base_tb_min = fill(Float32(150.0), full_channel_count)
        base_tb_max = fill(Float32(350.0), full_channel_count)
        required_valid_fraction = 0.5
    end

    if length(base_tb_min) < full_channel_count
        base_tb_min = vcat(base_tb_min, fill(Float32(150.0), full_channel_count - length(base_tb_min)))
    end
    if length(base_tb_max) < full_channel_count
        base_tb_max = vcat(base_tb_max, fill(Float32(350.0), full_channel_count - length(base_tb_max)))
    end

    selected_errors = base_errors[channel_indices]
    selected_use_flags = collect(use_flags[channel_indices])
    selected_tb_min = base_tb_min[channel_indices]
    selected_tb_max = base_tb_max[channel_indices]
    min_valid_fraction = clamp(required_valid_fraction, 0.1f0, 1.0f0)
    min_valid_channels = max(1, ceil(Int, min_valid_fraction * nchans))
    obs_error = repeat(selected_errors, 1, nobs)

    # Fill arrays from observations
    valid_obs_count = 0
    filter_reasons = Dict("missing_location" => 0, "invalid_location" => 0,
                         "zenith_too_large" => 0, "channel_mismatch" => 0,
                         "invalid_tb" => 0, "valid" => 0)

    for (i, obs) in enumerate(observations)
        # Extract geolocation - must be present
        if !haskey(obs, "latitude") || !haskey(obs, "longitude")
            filter_reasons["missing_location"] += 1
            continue  # Missing location data
        end

        lat_val = get(obs, "latitude", nothing)
        lon_val = get(obs, "longitude", nothing)

        if lat_val === nothing || lon_val === nothing
            filter_reasons["missing_location"] += 1
            continue
        end

        lat = Float32(lat_val)
        lon = Float32(lon_val)

        # Extract angles
        zenith = Float32(get(obs, "zenith_angle", 0.0))
        azimuth = Float32(get(obs, "azimuth_angle", 0.0))
        sol_zenith = Float32(get(obs, "solar_zenith", 0.0))
        sol_azimuth = Float32(get(obs, "solar_azimuth", 0.0))

        # Apply quality control filters
        if abs(lat) > 90.0 || abs(lon) > 180.0
            filter_reasons["invalid_location"] += 1
            continue  # Invalid location
        end

        if zenith > config.angle_limit
            filter_reasons["zenith_too_large"] += 1
            continue  # Zenith angle too large
        end

        # Extract brightness temperatures
        tb_values = get(obs, "brightness_temps", Float64[])
        if length(tb_values) != length(channels_list)
            filter_reasons["channel_mismatch"] += 1
            continue  # Mismatch in channel count
        end

        # Apply channel subset if needed
        if config.channel_subset !== nothing
            tb_values = tb_values[channel_indices]
        end

        # Convert to Float64 array if needed (Python objects)
        tb_array = try
            Float64.(tb_values)
        catch
            filter_reasons["invalid_tb"] += 1
            continue
        end

        # Basic physical QC - allow some missing/invalid channels
        # Set invalid channels to NaN for later filtering
        valid_mask = trues(nchans)
        for ch_idx in 1:nchans
            val = tb_array[ch_idx]
            if !isfinite(val)
                valid_mask[ch_idx] = false
            elseif val < selected_tb_min[ch_idx] || val > selected_tb_max[ch_idx]
                valid_mask[ch_idx] = false
            end
        end
        n_valid_channels = sum(valid_mask)

        # Require at least 50% of channels to be valid
        if n_valid_channels < min_valid_channels
            if filter_reasons["invalid_tb"] == 0 && config.verbose
                tb_min, tb_max = extrema(tb_array)
                @warn "Too few valid TB channels" n_valid=n_valid_channels n_total=nchans required=min_valid_channels tb_range=(tb_min, tb_max)
            end
            filter_reasons["invalid_tb"] += 1
            continue
        end

        # Set invalid channels to NaN (will be filtered in assimilation)
        tb_array[.!valid_mask] .= NaN

        filter_reasons["valid"] += 1
        valid_obs_count += 1

        # Store data
        brightness_temp[:, valid_obs_count] = convert(Vector{Float32}, tb_array)
        latitude[valid_obs_count] = lat
        longitude[valid_obs_count] = lon
        zenith_angle[valid_obs_count] = zenith
        azimuth_angle[valid_obs_count] = azimuth
        solar_zenith[valid_obs_count] = sol_zenith
        solar_azimuth[valid_obs_count] = sol_azimuth

        # Convert time
        base_time = get(obs, "base_time", (2000, 1, 1, 0, 0, 0))
        year = base_time[1] < 100 ? 2000 + base_time[1] : base_time[1]
        obs_time[valid_obs_count] = DateTime(year, base_time[2], base_time[3],
                                             base_time[4], base_time[5], base_time[6])
    end

    # Trim arrays to valid observations
    if valid_obs_count < nobs
        config.verbose && @info "Filtered to $valid_obs_count valid observations (from $nobs)" filter_reasons
        if valid_obs_count == 0
            @warn "All observations filtered out!" filter_reasons
        end
        brightness_temp = brightness_temp[:, 1:valid_obs_count]
        latitude = latitude[1:valid_obs_count]
        longitude = longitude[1:valid_obs_count]
        zenith_angle = zenith_angle[1:valid_obs_count]
        azimuth_angle = azimuth_angle[1:valid_obs_count]
        solar_zenith = solar_zenith[1:valid_obs_count]
        solar_azimuth = solar_azimuth[1:valid_obs_count]
        qc_flags = qc_flags[:, 1:valid_obs_count]
        obs_error = obs_error[:, 1:valid_obs_count]
        obs_time = obs_time[1:valid_obs_count]
    end

    # Create metadata
    metadata = Dict{String,Any}(
        "source_file" => filename,
        "format" => "Radiance BUFR",
        "satellite_id_code" => satellite_id_code,
        "instrument_id_code" => instrument_id_code,
        "total_observations" => valid_obs_count,
        "n_channels" => nchans,
        "channel_subset_applied" => config.channel_subset !== nothing,
        "channel_use_flags" => selected_use_flags,
        "channel_error_std" => selected_errors,
        "channel_tb_min" => selected_tb_min,
        "channel_tb_max" => selected_tb_max,
        "value_kind" => "brightnessTemperature",
        "source_value_kind" => source_value_kind,
        "qc_filters" => Dict(
            "max_zenith_angle" => config.angle_limit,
            "qc_threshold" => config.qc_threshold,
            "min_valid_fraction" => Float64(min_valid_fraction),
            "min_valid_channels" => min_valid_channels
        )
    )
    if settings !== nothing
        metadata["instrument_description"] = settings.description
    end

    return RadianceObservation(
        satellite_name,
        instrument_name,
        channels,
        brightness_temp,
        latitude,
        longitude,
        zenith_angle,
        azimuth_angle,
        solar_zenith,
        solar_azimuth,
        qc_flags,
        obs_error,
        obs_time,
        metadata
    )
end

"""
    convert_to_observation_data(radiance_obs::RadianceObservation) -> ObservationData

Convert RadianceObservation to generic ObservationData format for use in GSI core.

This flattens the multi-channel structure into individual observations.
"""
function convert_to_observation_data(radiance_obs::RadianceObservation)
    nchans, nobs = size(radiance_obs.brightness_temp)
    total_obs = nchans * nobs

    # Flatten arrays
    obs_values = vec(radiance_obs.brightness_temp)
    obs_errors = vec(radiance_obs.obs_error)
    qc_flags = vec(radiance_obs.qc_flags)

    # Expand location arrays to match flattened observations
    lons = repeat(radiance_obs.longitude, inner=nchans)
    lats = repeat(radiance_obs.latitude, inner=nchans)

    # Use channel number as "vertical" coordinate
    channels = repeat(radiance_obs.channels, outer=nobs)

    obs_locations = hcat(lons, lats, Float64.(channels))

    # Observation type code (radiance observations)
    # Use 200 + channel number to distinguish from conventional obs
    obs_types = fill(200, total_obs) .+ Int.(channels)

    # Expand time array
    obs_times = repeat(radiance_obs.obs_time, inner=nchans)

    # Metadata
    metadata = copy(radiance_obs.metadata)
    metadata["observation_category"] = "satellite_radiance"
    metadata["instrument"] = radiance_obs.instrument
    metadata["satellite"] = radiance_obs.satellite_id

    return ObservationData(obs_values, obs_errors, obs_locations, obs_types,
                          obs_times, qc_flags; metadata=metadata)
end

"""
    summarize_radiance_observations(radiance_obs::RadianceObservation) -> Dict

Generate summary statistics for radiance observations.
"""
function summarize_radiance_observations(radiance_obs::RadianceObservation)
    nchans, nobs = size(radiance_obs.brightness_temp)

    # Retrieve per-channel metadata when available
    channel_use_flags = get(radiance_obs.metadata, "channel_use_flags", fill(Int8(1), nchans))
    channel_error_std = get(radiance_obs.metadata, "channel_error_std", fill(Float32(2.0), nchans))

    # Compute statistics per channel
    channel_stats = Dict{Int,Dict{String,Float64}}()
    for (i, ch) in enumerate(radiance_obs.channels)
        tb_chan = radiance_obs.brightness_temp[i, :]
        valid_mask = .!isnan.(tb_chan)
        valid_vals = tb_chan[valid_mask]
        valid_count = count(valid_mask)
        idx = i <= length(channel_error_std) ? i : length(channel_error_std)
        uf_idx = i <= length(channel_use_flags) ? i : length(channel_use_flags)
        error_std = Float64(channel_error_std[idx])
        use_flag = Float64(channel_use_flags[uf_idx])

        if isempty(valid_vals)
            channel_stats[ch] = Dict(
                "mean_tb" => NaN,
                "std_tb" => NaN,
                "min_tb" => NaN,
                "max_tb" => NaN,
                "obs_error_std" => error_std,
                "obs_error_var" => error_std^2,
                "use_flag" => use_flag,
                "valid_count" => 0.0,
                "valid_fraction" => 0.0
            )
        else
            stats = Dict{String,Float64}(
                "mean_tb" => mean(valid_vals),
                "std_tb" => std(valid_vals),
                "min_tb" => minimum(valid_vals),
                "max_tb" => maximum(valid_vals),
                "obs_error_std" => error_std,
                "obs_error_var" => error_std^2,
                "use_flag" => use_flag,
                "valid_count" => Float64(valid_count),
                "valid_fraction" => valid_count / nobs
            )
            channel_stats[ch] = stats
        end
    end

    summary = Dict{String,Any}(
        "satellite" => radiance_obs.satellite_id,
        "instrument" => radiance_obs.instrument,
        "n_observations" => nobs,
        "n_channels" => nchans,
        "channels" => radiance_obs.channels,
        "channel_statistics" => channel_stats,
        "spatial_coverage" => Dict(
            "lat_range" => (minimum(radiance_obs.latitude), maximum(radiance_obs.latitude)),
            "lon_range" => (minimum(radiance_obs.longitude), maximum(radiance_obs.longitude))
        ),
        "geometry" => Dict(
            "zenith_angle_range" => (minimum(radiance_obs.zenith_angle), maximum(radiance_obs.zenith_angle)),
            "mean_zenith_angle" => mean(radiance_obs.zenith_angle)
        ),
        "time_range" => (minimum(radiance_obs.obs_time), maximum(radiance_obs.obs_time))
    )

    return summary
end

end # module RadianceBUFRReader
